En general, cuantificar incertidumbre es computacionalmente costoso y/o metodológicamente más desafiante que no hacerlo. Esto lleva a ignorar la incertidumbre cuando tenemos modelos muy complejos, como deep neural networks, o sets de datos muy grandes. Pero además, cuantificar o no incertidumbre en la estimación también depende del contexto. Cuando simplemente necesitamos un modelo que prediga bien y no nos preocupa mucho cómo lo hace, o cuando es difícil juzgar la performance de un modelo, cuantificar incertidumbre suele ser secundario. Por ejemplo, ¿cómo podríamos evaluar la incertidumbre de una respuesta de un modelo de lenguaje?; ¿cuánto nos preocupa la incertidumbre en la predicción de un modelo que genera sugerencias de contenido en Netflix?
Optimización
Minimización de funciones de costo
El enfoque más sencillo es minimizar una medida de discrepancia entre los datos y la predicción.
En general, llamamos a esta medida función de costo (loss function).
Consideremos un modelo muy simple: \[\hat{y} = \lambda\]
Con una sola observación, \(y_1\), la discrepancia se minimiza en \[\lambda = y_1\]
Tomamos como función de costo el error absoluto: \(|y_i – \hat{y_i}|.\)
Code
lambda_seq <-seq(0, 30, by =0.25)# Error absolutoy1 <-11d1 <-abs(y1 - lambda_seq)plot(d1 ~ lambda_seq, type ="l", col =2, xlab =expression(lambda), ylab ="Error absoluto")text(2.5, 12, "y = 11", col =2)y2 <-5d2 <-abs(y2 - lambda_seq)lines(d2 ~ lambda_seq, col =3)text(2.5, 14, "y = 5", col =3)y3 <-20d3 <-abs(y3 - lambda_seq)lines(d3 ~ lambda_seq, col =4)text(2.5, 16, "y = 20", col =4)
Con más de una observación, requerimos una métrica global.
Error Absoluto Medio (MAE): \[\frac{1}{N} \sum_{i=1}^N |y_i – \hat{y_i}|\]
Code
plot(d1 ~ lambda_seq, type ="l", col =2, xlab =expression(lambda), ylab ="Error absoluto")text(2.5, 12, "y = 11", col =2)lines(d2 ~ lambda_seq, col =3)text(2.5, 14, "y = 5", col =3)lines(d3 ~ lambda_seq, col =4)text(2.5, 16, "y = 20", col =4)# MAEmae <-rowMeans(cbind(d1, d2, d3))lines(mae ~ lambda_seq, lwd =2, lty =2)text(12, 14, "MAE")abline(v =median(c(y1, y2, y3)), lty =2)
¿Promedio o suma?
El promedio y la suma son operaciones equivalentes, ya que generalmente lo relevante es la forma de la función de costo, no su valor absoluto. El promedio es una transformación lineal de la suma (se multiplica por 1/N), y ninguna transformación lineal altera la forma. Es decir, podemos sumar los errores y multiplicarles cualquier constante o sumarles cualquier constante y no tendremos problemas.
El MAE se minimiza en la mediana.
No es derivable (problemático).
e1 <- (y1 - lambda_seq) ^2plot(e1 ~ lambda_seq, type ="l", col =2, xlab =expression(lambda), ylab ="Error cuadrático")text(2.5, 150, "y = 11", col =2)e2 <- (y2 - lambda_seq) ^2lines(e2 ~ lambda_seq, col =3)text(2.5, 180, "y = 5", col =3)e3 <- (y3 - lambda_seq) ^2lines(e3 ~ lambda_seq, col =4)text(2.5, 210, "y = 20", col =4)# error absoluto mediamse <-rowMeans(cbind(e1, e2, e3))lines(mse ~ lambda_seq, lwd =2, lty =2)text(12, 180, "MSE")abline(v =mean(c(y1, y2, y3)), lty =2)
Code
# Y da lo mismo si sumamos los errores en vez de promediar:tse <-rowSums(cbind(e1, e2, e3))plot(tse ~ lambda_seq, type ="l", lwd =2, lty =2, xlab =expression(lambda), ylab ="Error cuadrático (suma)")abline(v =mean(c(y1, y2, y3)), lty =2)
Detalles
El error cuadrático se minimiza en la media de los datos, y, al igual que con el MAE, da lo mismo si sumamos o promediamos los errores individuales.
Una función de costo elegante: la función de verosimilitud
Los modelos estadísticos definen de manera explícita la relación entre las predicciones y los datos.
Los datos son considerados realizaciones de una variable aleatoria, cuya distribución se modela. Ejemplo:
lambda_seq <-seq(0, 30, by =0.1)l1 <-dpois(y1, lambda_seq)plot(l1 ~ lambda_seq, type ="l", col =2, ylim =c(0, 0.2),ylab =bquote("Pr(y | "* lambda *")"), xlab =expression(lambda),main ="Función de verosimilitud")text(26, 0.14, expression(y[1] ==11), col =2)l2 <-dpois(y2, lambda_seq)lines(l2 ~ lambda_seq, col =3)text(26, 0.12, expression(y[2] ==5), col =3)l3 <-dpois(y3, lambda_seq)lines(l3 ~ lambda_seq, col =4)text(26, 0.10, expression(y[3] ==20), col =4)# verosimilitud conjunta:lfull <- l1 * l2 * l3lines(lfull ~ lambda_seq, lwd =2, lty =2)text(26, 0.18, "Verosimilitud\nconjunta")
¡Es plana!
La escalamos
Code
lambda_seq <-seq(0, 30, by =0.1)l1 <-dpois(y1, lambda_seq)plot(l1 ~ lambda_seq, type ="l", col =2, ylim =c(0, 0.2),ylab =bquote("Pr(y | "* lambda *")"), xlab =expression(lambda),main ="Función de verosimilitud")text(26, 0.14, expression(y[1] ==11), col =2)l2 <-dpois(y2, lambda_seq)lines(l2 ~ lambda_seq, col =3)text(26, 0.12, expression(y[2] ==5), col =3)l3 <-dpois(y3, lambda_seq)lines(l3 ~ lambda_seq, col =4)text(26, 0.10, expression(y[3] ==20), col =4)# verosimilitud conjunta:lfull <- l1 * l2 * l3a <-diff(lambda_seq)[1]lfull_norm <-normalize_dens(lfull, a) # sólo la escalamos, para que se vealines(lfull_norm ~ lambda_seq, lwd =2, lty =2)text(26, 0.18, "Verosimilitud\nconjunta")
La función de verosimilitud se va volviendo más puntuda a medida que agregamos datos:
Code
y <-c(y1, y2, y3, 4, 10, 12) # lo extendemos# vector con verosimilitud para cada valor de lambda_seql <-numeric(length(lambda_seq))for (i in1:length(lambda_seq)) { l[i] <-prod(dpois(y, lambda_seq[i]))}l_norm <-normalize_dens(l, diff(lambda_seq)[1])plot(l_norm ~ lambda_seq, type ="l",xlab =expression(lambda), ylab =expression(Pr(bold(y) ~"|"~ lambda)))lines(lfull_norm ~ lambda_seq, col =2)text(26, 0.25, "N = 6")text(26, 0.22, "N = 3", col =2)
Optimización:
Buscamos el valor de \(\lambda\) que maximice la verosimilitud: estimador de máxima verosimilitud.
(maximum likelihood estimate, MLE)
Generalmente se trabaja con la log-verosimilitud.
Recordemos que
\[\log(a \times b) = \log(a) + \log(b).\]
Entonces, la log-verosimilitud conjunta se define así:
ll3 <-log(lfull) # la que calculamos con 3 datosll6 <-log(l) # la que calculamos con 6 datosrr <-c(ll3, ll6)rr <-range(rr[is.finite(rr)])plot(ll6 ~ lambda_seq, type ="l",xlab =expression(lambda), ylab =expression(log(Pr(bold(y) ~"|"~ lambda))),ylim = rr)lines(ll3 ~ lambda_seq, col =2)text(26, -80, "N = 6")text(26, -40, "N = 3", col =2)
Y si la negamos, nos interesa minimizarla:
Code
plot(-ll6 ~ lambda_seq, type ="l",xlab =expression(lambda), ylab =expression(-log(Pr(bold(y) ~"|"~ lambda))),ylim = rr[2:1] * (-1))lines(-ll3 ~ lambda_seq, col =2)text(26, 80, "N = 6")text(26, 40, "N = 3", col =2)
En un contexto de optimización, los modelos estadísticos son una forma elegante de elegir la función de costo, definida como la -log verosimilitud.
Máxima verosimilitud en 2D (ejemplo)
Para tener un ejemplo real, podemos modelar el número de incendios por verano en función de un índice de peligrosidad de incendios, llamado FWI (datos de Barberá et al. 2025).
\[
\begin{aligned}
y_i &\sim \text{poisson}(\lambda_i) \\
\lambda_i &= \text{exp}(\alpha + \beta \ FWI_i) \\
\\
y_i&: \text{número de incendios} \ge 10 \text{ha en el verano} \ i \\
FWI_i&: \text{Fire Weather Index en el verano} \ i
\end{aligned}
\]
\(\lambda\) no se estima, es una cantidad derivada.
Sólo necesitamos estimar \(\alpha\) y \(\beta\).
Cargamos datos y graficamos.
Code
datos <-read.csv(here::here("R_curso_modelos", "datos", "barbera_data_fire_total_climate.csv"))ggplot(datos, aes(x = fwi, y = fires)) +geom_point(color =rgb(0, 0, 0, 0.8), size =2) +ylab("y") +xlab("FWI") +nice_theme()
Función de verosimilitud, ahora depende de \(\alpha\) y \(\beta\)
(uso interactivo, vaya a R).
Code
like_fire <-function(alpha, beta, log = T) { # "like" por likelihood (sería raro decirle "vero");# "fire" porque es la función para estos datos de fuego, no cualquier # verosimilitud.# Calculamos la media, lambda: lambda <-exp(alpha + beta * datos$fwi)# usando lambda evaluamos la verosimilitud de cada observación, # en escala log: like_pointwise <-dpois(datos$fires, lambda, log = T)# la log-verosimilitud conjunta es la suma de las log-verosimilitudes por # observación: like <-sum(like_pointwise)if (log) return(like)elsereturn(exp(like))}
Ahora podemos darle valores a \(\alpha\) y \(\beta\) y evaluar la verosimilitud.
También graficamos la curva definida por esos valores.
(Uso interactivo, vaya a R.)
Code
## Podemos elegir valores a mano# alpha <- log(2)# beta <- 0.01# O más interesante, muestrear en cierto rangoalpha <-runif(1, log(0.01), log(30))beta <-runif(1, 0, 0.2)like <-like_fire(alpha, beta)# Label para el gráficolabel_text <-as.character(as.expression(bquote(L(y ~"|"~ alpha == .(round(alpha, 3)) ~","~ beta == .(round(beta, 3))) == .(like, 3)) ))# Plot de puntos + curvaggplot(datos, aes(x = fwi, y = fires)) +geom_point(color =rgb(0, 0, 0, 0.8), size =3) +stat_function(fun =function(x) exp(alpha + beta * x), color ="blue", linewidth =1) +ylab("y") +xlab("FWI") +annotate("text", x =12, y =25,label = label_text,parse =TRUE, hjust =0.5) +ylim(min(datos$fires -1), max(datos$fires +1)) +xlim(min(datos$fwi), max(datos$fwi)) +nice_theme()
Mejor que probar valores al azar es buscar de forma más ordenada. El método más simple, aunque ineficiente, es buscar en una grilla.
Code
# Creamos grilla de valores de alpha y beta.side <-100# cantidad de valores por ladogrilla <-expand.grid(alpha =seq(log(0.01), log(30), length.out = side),beta =seq(0, 0.3, length.out = side))# Evaluamos la verosimilitud para cada valorgrilla$loglike <-NAsize <- side ^2# cantidad de valores en la grillafor (i in1:size) { grilla$loglike[i] <-like_fire(grilla$alpha[i], grilla$beta[i]) }grilla$like <-exp(grilla$loglike)
Visualizamos la grilla y evaluamos puntos.
Code
# Elegimos un punto para evaluarrow <-sample(1:size, 1)alpha <- grilla$alpha[row]beta <- grilla$beta[row]# Grillap1 <-ggplot(grilla, aes(x = alpha, y = beta, fill = like)) +geom_tile() +geom_point(data = grilla[row, , drop = F]) +scale_fill_viridis_c(option ="viridis", name ="Verosimilitud",label =scientific_format(digits =2)) +labs(x =expression(alpha), y =expression(beta)) +theme_minimal(base_size =14)# Plot de puntos + curvap2 <-ggplot(datos, aes(x = fwi, y = fires)) +geom_point(color =rgb(0, 0, 0, 0.8), size =2) +stat_function(fun =function(x) exp(alpha + beta * x), color ="blue", linewidth =1) +ylab("y") +xlab("FWI") +nice_theme()# Unimos(p1 | p2)
Y como estimador puntual tenemos el valor con mayor verosimilitud:
La diferencia entre los enfoques puede parecer enorme, y ha habido mucho debate sobre cuál es el menos feo. Sin embargo, desde una perspectiva práctica, los enfoques tienen mucho en común, salvo quizás cuando se trata de la selección de modelos. En particular, si se aplican correctamente, suelen producir resultados que difieren menos entre sí de lo que probablemente difieren los modelos analizados con respecto a la realidad.
[Wood, Simon. 2015. Core Statistics. Cambridge University Press.]
Inferencia bayesiana
Nuestro conocimiento sobre los parámetros se representa con distribuciones de probabilidad. Es decir, los tratamos como variables aleatorias.
Para pensar en un modelo y la distribución de sus parámetros no necesitamos estrictamente tener datos.
Si tenemos información nueva (datos), actualizamos nuestro conocimiento. En este caso, llamamos previa a la distribución que antecede a la evidencia, posterior a la distribución actualizada, y la evidencia nueva está codificada en la verosimilitud.
Teorema de Bayes:
\[\Pr(A \mid B) = \frac{\Pr(B \mid A) \times \Pr(A)} {\Pr(B)}\]
# Creamos grilla de valores de alpha y beta.side <-100# cantidad de valores por ladoalpha_seq <-seq(log(0.01), log(30), length.out = side)beta_seq <-seq(-0.5, 0.5, length.out = side)grilla <-expand.grid(alpha = alpha_seq,beta = beta_seq,loglike =NA, # para llenar luegolike =NA,prior =NA,post =NA)# Evaluamos la verosimilitud para cada valorsize <- side ^2# cantidad de valores en la grillafor (i in1:size) { grilla$loglike[i] <-like_fire(grilla$alpha[i], grilla$beta[i]) }grilla$like <-exp(grilla$loglike)# Definimos y evaluamos la previamu_a <-0; sigma_a <-1mu_b <-0; sigma_b <-0.1for (i in1:size) { grilla$prior[i] <-dnorm(grilla$alpha[i], mu_a, sigma_a) *dnorm(grilla$beta[i], mu_b, sigma_b)}# Posterior grilla$post <- grilla$like * grilla$prior# Normalizamos para poder visualizarg <- grillaa <-diff(alpha_seq)[1] *diff(beta_seq)[1] # para normalizarfor (v inc("like", "prior", "post")) { g[, v] <-normalize_dens(grilla[, v], a)}# Graficamos por separado, para apreciar las distintas escalasp1 <-ggplot(g, aes(x = alpha, y = beta, fill = like)) +geom_tile() +scale_fill_viridis() +theme(legend.title =element_blank()) +xlab(expression(alpha)) +ylab(expression(beta)) +ggtitle("Verosimilitud") +nice_theme()p2 <-ggplot(g, aes(x = alpha, y = beta, fill = prior)) +geom_tile() +scale_fill_viridis() +theme(legend.title =element_blank()) +xlab(expression(alpha)) +ylab(expression(beta)) +ggtitle("Previa") +nice_theme()p3 <-ggplot(g, aes(x = alpha, y = beta, fill = post)) +geom_tile() +scale_fill_viridis() +theme(legend.title =element_blank()) +xlab(expression(alpha)) +ylab(expression(beta)) +ggtitle("Posterior") +nice_theme()(p1 | p2 | p3)# + plot_layout(nrow = 2)
¿Cómo definir en R una función que evalúe la densidad posterior?
Code
post_fire <-function(alpha, beta, mu_a =0, sigma_a =1, # parámetros de la previamu_b =0, sigma_b =1, log = T) { ## Log Verosimilitud (mismo código que en like_fire)# Calculamos la media, lambda: lambda <-exp(alpha + beta * datos$fwi)# usando lambda evaluamos la verosimilitud de cada observación, # en escala log: like_pointwise <-dpois(datos$fires, lambda, log = T)# la log-verosimilitud conjunta es la suma de las log-verosimilitudes por # observación: like <-sum(like_pointwise)## Log Previas lprior_a <-dnorm(alpha, mean = mu_a, sd = sigma_a, log = T) lprior_b <-dnorm(beta, mean = mu_b, sd = sigma_b, log = T)## Log Posterior lpost <- like + lprior_a + lprior_bif (log) return(lpost)elsereturn(exp(lpost))}
Trabajo Práctico 5.2
Actividades 2.
Ajuste frecuentista hegemónico a mano, usando optim
La práctica más usual para calcular valores P o intervalos de confianza en el enfoque frecuentista implica suponer que la función de log verosimilitud tiene forma cuadrática, centrada en el máximo. Esto se cumple si los parámetros no se encuentran cerca de un límite (e.g., una varianza muy cercana a cero) y si el N es grande en relación a la cantidad de parámetros.
La información sobre la curvatura de la log verosimilitud se encuentra en la matriz hessiana, que es la matriz de derivadas parciales de segundo orden de la log verosimilitud evaluadas en el máximo. optim la calcula si le indicamos hessian = TRUE:
Code
opt <-optim(par =c(0, 0.1), # vector inicial de parámetros, donde comienza la búsquedafn = like_fire_opt, # función a optimizarmethod ="BFGS",hessian =TRUE)opt
Si repitiéramos el estudio (incluyendo este análisis) infinitas veces, siempre tomando un N grande, obtendríamos una muestra de estimadores de máxima verosimilitud (una por estudio), los cuales seguirían una Normal Multivariada alrededor del valor verdadero de los parámetros. A partir de esto, se desprende un método para construir intervalos de confianza y calcular valores P, que utilizamos a continuación:
Con esta matriz, podemos construir una tabla de resumen análoga a las que arrojan funciones como glm.
Code
sds <-sqrt(diag(I)) # error estándar de los estimadoreszs <- opt$par / sds # valor zps <- (1-pnorm(abs(zs), mean =0, sd =1)) *2# valor P, a dos colassumm_table <-data.frame(par =c("alpha", "beta"),est = opt$par,sd = sds,z = zs,p = ps)names(summ_table) <-c("Parameter", "Estimate", "Std. Error", "z value", "Pr(>|z|)")summ_table
Parameter Estimate Std. Error z value Pr(>|z|)
1 alpha 0.08094103 0.29263313 0.2765956 7.820907e-01
2 beta 0.16529256 0.02000867 8.2610459 2.220446e-16
Code
# Comparamos con glmglm_fit <-glm(fires ~ fwi, data = datos, family ="poisson")summary(glm_fit)
Call:
glm(formula = fires ~ fwi, family = "poisson", data = datos)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.08042 0.29267 0.275 0.783
fwi 0.16532 0.02001 8.261 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 122.571 on 23 degrees of freedom
Residual deviance: 54.035 on 22 degrees of freedom
AIC: 150.78
Number of Fisher Scoring iterations: 4
Máxima verosimilitud penalizada: casi una previa
Si en vez de concentrarnos en obtener la distribución posterior completa nos conformamos únicamente con obtener el modo, como estimador puntual, podríamos tomar la -log densidad posterior y minimizarla. La -log densidad posterior, es la suma de la -log verosimilitud, la -log previa y el -log de la constante normalizadora (\(C\)):
En donde ignoramos \(C\) porque no varía en función de los parámetros, y el sufijo \(MAP\) viene de maximum a posteriori, el modo de la posterior.
Esta expresión es muy similar a un estimador de máxima verosimilitud penalizada (mínima -log verosimilitud), en donde se introduce un término de penalidad \(K\):
Con \(PML\) por penalized maximum likelihood.
La penalidad \(K\) suele utilizarse para evitar que un modelo muy complejo se sobreajuste a los datos. En un contexto de regresión múltiple, con muchas predictoras, \(K\) se define como una función creciente del valor absoluto de los coeficientes (\(\beta_j\)). Por ejemplo, en LASSO (Least Absolute Shrinkage and Selection Operator), \(K\) se define así:
\[K_{LASSO} = \lambda \ \sum_{j=1}^J |\beta_j|,\]
siendo \(\lambda > 0\) un parámetro de penalización usualmente escogido mediante validación cruzada. Esta forma de penalización también es conocida como L1. Desde una perspectiva bayesiana, introducir esta penalidad equivale a definir una previa de Laplace (doble exponencial) centrada en cero sobre los coeficientes:
De todos modos, recordemos que estimar simplemente el modo de la posterior no es una práctica usual en el enfoque bayesiano. Además de esto, la diferencia entre el estimador puntual bayesiano y la verosimilitud penalizada que en el primero la amplitud de la previa, definida por \(\lambda\), se escoge sin mirar los datos. En cambio, en un contexto de optimización, se estima utilizando los datos (por ejemplo, mediante validación cruzada). Pero también hay bayesianos que escogen el valor de \(\lambda\) en base a los datos, y este enfoque se llama Bayesianismo empírico.
Otra forma de penalización es Ridge Regression, también llamada L2:
Estas formas de penalización son muy utilizadas en machine learning para ajustar modelos de muchos parámetros previniendo el sobreajuste. También se suele usar la elastic net, una combinación lineal de Ridge y LASSO:
\[K_{EN} = \lambda \ [\alpha \sum_{j=1}^J \beta_j^2 + (1-\alpha) \sum_{j=1}^J |\beta_j|],\] con \(\alpha \in (0, 1)\). Esto equivale a una previa definida como una mezcla entre una distribución de Laplace y una Normal, en donde \(\alpha\) regula el peso de cada distribución.
Y otra forma de penalización es la Bridge:
\[K_{Bridge} = \lambda \sum_{j=1}^J |\beta_j|^p,\] que se reduce a LASSO cuando \(p = 1\) y a Ridge con \(p = 2\).
La estrategia de agregar términos de penalización a la -log verosimilitud para evitar el sobreajuste también es conocida como regularización.